NPS-MA-95-003 


NAVAL  POSTGRADUATE  SCHOOL 
Monterey,  California 


AN  ALGORITHM  FOR  COMPUTING  THE 
STATIONARY  DISTRIBUTION  OF  A  DISCRETE-TIME 
BIRTH-AND-DEATH  PROCESS  WITH  BANDED 
INFINITESIMAL  GENERATOR 

by 

Carlos  F.  Borges 
Craig  S.  Peters 

April  1995 

Report  for  Period:  November  1994  -  January  1995 


Approved  for  public  release;  distribution  is  unlimited 

Prepared  for:  Naval  Postgraduate  School 
Monterey,  CA  93943 


dtxc  QSMm 


19950703  054 


NAVAL  POSTGRADUATE  SCHOOL 
MONTEREY,  CA  93943 


Rear  Admiral  T.A.  Mercer 
Superintendent 


Harrison  Shull 
Provost 


This  report  was  prepared  in  conjunction  with  research  conducted  for  the  Naval  Postgraduate  School  and 
funded  by  the  Naval  Postgraduate  School. 

Reproduction  of  all  or  part  of  this  report  is  authorized. 

This  report  was  prepared  by: 


&  & 


REPORT  DOCUMENTATION  PAGE 


Form  Approved 
OM8  No.  0704-0188 


PuoJic  reoortinq  burden  for  this  collection  of  information  is  estimated  to  average  l  hour  per  response,  including  the  time  for  reviewing  instructions,  searching  existing  data  sources, 
aarherina  and  maintaining  the  data  needed,  and  completing  and  reviewing  the  collection  of  information.  Send  comments  regarding  this  burden  estimate  or  any  other  aspect  of  this 
collection  of  information,  including  suggestions  for  reducing  this  burden,  to  Washington  Headquarters  Services,  Directorate  tor ,n f 9°t.r*l!ons  and  5  Jetferv:,n 

Davis  Highway.  Suite  1204.  Arlington,  VA  22202-4302.  and  to  the  Office  of  Management  and  Budget,  Paperwork  Reduction  Project  (Q7Q4-Q 188).  Washington.  DC  20503. 


1.  AGENCY  USE  ONLY  (Leive  blank )  2.  REPORT  DATE 

April  8,  1995 


3.  REPORT  TYPE  AND  DATES  COVERED 

Technical  Report  Nov  1994  -  Jan  1995 


4.  TITLE  AND  SUBTITLE 

An  Algorithm  for  Computing  the  Stationary  Distribution  of  a  Discrete- 
Time  Birth-and-Death  Process  with  Banded  Infinitesimal  Generator 


6.  AUTHOR(S) 


Carlos  F.  Borges  and  Craig  S.  Peters 


7.  PERFORMING  ORGANIZATION  NAME(S)  AND  AOORESS(ES) 

Naval  Postgraduate  School 
Monterey,  CA  93943-5000 


8.  PERFORMING  ORGANIZATION 
REPORT  NUMBER 

NPS-MA-95-003 


9.  SPONSORING /MONITORING  AGENCY  NAME(S)  AND  ADDRESS(ES) 


10.  SPONSORING /MONITORING 
AGENCY  REPORT  NUMBER 


Naval  Postgraduate  School 
Monterey,  CA  93943 


i*hJuviewsE expressed  in  this  report  are  those  of  the  authors  and  do  not  reflect  the  official  policy  or 
position  of  the  Department  of  Defense  or  the  United  States  Government. 


12a.  DISTRIBUTION /AVAILABILITY  STATEMENT 

Approved  for  public  release;  distribution  is  unlimited. 


12b.  DISTRIBUTION  CODE 


13.  ABSTRACT  (Maximum  200  words) 

We  develop  an  algorithm  for  computing  approximations  to  the  stationary  distribution  of  a  discrete¬ 
time  birth-and-death  process  provided  that  the  infinitesimal  generator  is  a  banded  matrix.  We  begin 
by  computing  stationary  distributions  for  processes  whose  infinitesimal  generators  are  Hessenburg. 
Our  derivation  in  this  special  case  is  different  than  the  classical  one  but  leads  to  the  same  result. 

We  then  show  how  to  extend  these  ideas  to  get  approximations  when  the  infinitesimal  generator  is 
banded  (or  half-banded). 


14.  SUBJECT  TERMS  ,  ..... 

algorithm  for  computing  approximations  to  the  stationary  distribution 


15.  NUMBER  OF  PAGES 

14 


15.  PRICE  CODE 


17.  SECURITY  CLASSIFICATION  I  18.  SECURITY  CLASSIFICATION  19.  SECURITY  CLASSIFICATION  J  20.  LIMITATION  OF  ABSTRACT 

OF  REPORT  OF  THIS  PAGE  OF  ABSTRACT  I 


[UNCLASSIFIED 

NSN  7540-01-280-5500 


UNCLASSIFIED 


UNCLASSIFIED 


Standard  Form  298  (Rev  2-89) 
Prescribed  by  ANSI  Svb  239- '8 
298-102  _ 


An  Algorithm  for  Computing  the  Stationary 
Distribution  of  a  Discrete  Birth-and-Death  Process 
with  Banded  Infinitesimal  Generator 

Technical  Report:  NPS-MA-95-003 

Carlos  F.  Borges  *  Craig  S.  Peters 
January  5,  1995 


Abstract 

We  develop  an  algorithm  for  computing  approximations  to  the  sta¬ 
tionary  distribution  of  a  discrete  birth-and-death  process  provided  that 
the  infinitesimal  generator  is  a  banded  matrix.  We  begin  by  comput¬ 
ing  stationary  distributions  for  processes  whose  infinitesimal  generators 
are  Hessenberg.  Our  derivation  in  this  special  case  is  different  than  the 
classical  one  but  leads  to  the  same  result.  We  then  show  how  to  extend 
these  ideas  to  get  approximations  when  the  infinitesimal  generator  is 
banded  (or  half-banded). 


1  Introduction 

A  birth  and  death  process  is  characterized  by  a  population  of  individuals 
whose  number  changes  according  to  the  outcome  of  two  other  stochastic  pro¬ 
cesses  consisting  of  births  which  increase  the  population  and  deaths  which 
decrease  the  population.  The  transition  probabilities  for  these  two  processes 
can,  in  general,  depend  on  both  time  and  population  size.  The  model  for 
this  stochastic  dynamical  system  is  usually  described  by  a  master  equation 
for  the  transition  probability  for  population  size  at  a  given  time. 

Let  N(t )  denote  the  population  size  at  time  t  and  define  the  transition 
probability  for  N(t)  as  P(n,t)  =  Pr{iV(t)  =  n|W(F)  =  n'}.  The  transi¬ 
tion  probabilities  for  births  and  deaths  are  usually  modeled  as  Markovian 
processes  by  assuming  that 
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r{n,t )  =  Pr{N(t  +  dt)  =  k  +  l\N(t)  =  k} 
l(n,t )  =  Pr{N(t  +  dt)  =  k  -  ljN(t)  =  k} 

where  dt  represents  the  fundamental  time  unit.  Time  can  be  treated  as  either 
a  discrete  or  continuous  variable.  In  many  situations  it  is  reasonable  and 
convenient  to  model  these  transition  probabilities  as  being  time  independent 
and  to  assume  that  the  probability  that  more  than  a  single  birth  or  death 
occurs  in  the  fundamental  time  unit  is  zero.  With  these  assumptions  the 
master  equation  for  P(n,  t )  can  be  written 


P(n,t  +  dt)  =  r(n  —  l)P(n  —  l,i)  + /(»  +  1)P(»+ l,t) + 

+  [l-(r(n)  +  Z(n))]P(»,t)  (1) 

The  first  approach  to  investigating  1  is  to  assume  that 


lim  P(n,  t  +  dt)  —  P(n ,  t)  =  0 

£— KX) 

and  write  the  equation  for  the  stationary  transition  probabilities,  p(n),  the 
probability  that  the  population  will  eventually  stabilize  at  n  individuals,  as 


0  =  — (r(n)  +  l(n))p(n)  +  r(n  -  1  )p(n  -  1)  +  l(n  +  1  )p{n  +  1)  (2) 

This  stochastic  balance  equation  leads  to  the  infinitesimal  generator  for 
the  discrete  time  Markov  process  which  has  the  following  stationary  distri¬ 
bution 


p(n)  =  p(0)  f[  r(^y--  (3) 

In  this  paper  we  develop  an  algorithm  for  computing  the  solution  to  this 
problem  and  show  how  it  can  be  generalized  to  compute  approximate  solu¬ 
tions  for  birth-and-death  processes  where  both  multiple  births  and  multiple 
deaths  can  occur  in  the  fundamental  time  unit.  These  methods  are  devel¬ 
oped  by  converting  stochastic  balance  equations  like  2  to  matrix  form  and 
applying  techniques  from  linear  algebra. 

The  virtue  of  the  linear  algebra  approach  is  twofold.  First,  the  use  of 
finite  precision  calculations  in  linear  algebra  is  well  understood  and  many 
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algorithms  have  been  developed  that  can  control  the  ill  effects  of  round¬ 
off  and  other  errors.  Second,  generalizations  of  the  simple  birth  and  death 
model  give  rise  to  generalizations  of  the  master  equations  1  and  2  for  which 
solutions  are  not  known.  In  this  case,  the  methods  we  develop  are  still 
applicable.  To  motivate  what  follows,  we  now  show  how  3  results  from 
solving  an  appropriate  linear  system. 


2  A  Matrix  Formulation 

We  are  interested  in  finding  the  stationary  distribution  p(n)  which  satisfies 
equation  2.  We  begin  by  converting  2  to  matrix  form.  In  what  follows  vectors 
will  be  denoted  by  lower-case  bold  Roman  letters  and  will  be  assumed  to  be 
column  vectors.  We  shall  denote  by  0  the  vector  of  all  zeros,  by  e  the  vector 
of  all  ones,  and  by  e;  the  i’th  axis  vector,  a  vector  whose  i’th  element  is  one 
and  all  others  are  zero.  The  sizes  of  these  vectors,  when  they  appear,  shall 
be  taken  from  context. 

To  proceed,  we  represent  the  stationary  distribution  p(n )  as  an  infinite 
column  vector  p  whose  i’th  element  p,  =  p(i  —  1)  where  i  =  1,2, ...,  oo.  Note 
that  if  there  were  a  maximum  population  size  N  then  p  would  be  an  element 
of  iftjV+1.  In  general,  however,  this  is  not  the  case. 

The  infinitesimal  generator  matrix  QT  is  an  infinite  tridiagonal  matrix 
with  entries 


— (r(i  -  1)  +  l(i  -  1))  if  i  =  j 

r(j-l)  if*-l=j  m 

l(j  -  1)  if  i  +  1  =  j 

0  otherwise 

With  these  definitions,  equation  2  becomes 

QT  P  =  0  (5) 

And  we  see  that  QT  must  be  rank- deficient  and  p  in  its  null-space  (or  equiv¬ 
alently,  p  is  an  eigenvector  of  QT  associated  with  the  eigenvalue  A  =  0). 
Equation  5  implies  that  any  element  of  the  null-space  of  0  ^  solves  the  equa¬ 
tion.  To  get  the  stationary  distribution  we  normalize  using  the  law  of  total 
probability,  so  that  eTp  =  1. 
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3  Truncated  Solutions 


It  is  not  generally  possible  to  solve  infinite  systems  of  equations,  or  find 
eigenvectors  of  infinite  matrices,  so  we  consider  truncating  the  infinite  sys¬ 
tem  of  equations.  Graphically,  we  truncate  the  infinite  system  by  by  parti¬ 
tioning  it  in  the  following  way 


o 

X 

X 

x  **.  *•.  *** 

O 

X 

X 

X 

X 

X 

X  X 

X 

p(0) 

K1) 


P(n~  !) 

p{n) 


=  0 


and  then  dropping  all  but  the  first  n  equations. 

Letting  Q ^  be  the  first  principal  n  x  n  sub-matrix  of  QT,  and  pn  be  the 
first  principal  n  vector  of  p,  we  get 


<2nPn  =  ~l(n)p(n)en 

In  effect,  this  is  a  matrix  representation  of  the  first  n  equations  from  5. 
Letting  pn  =  l(n)p(n)ta  and  rearranging  yields 

Qn  fn  =  -e„  (6) 

So,  provided  that  Q ^  is  non-singular  and  that  l(n)p(n)  /  0,  we  can  solve 
directly  for  a  scalar  multiple  of  the  truncated  stationary  distribution  pn. 

In  particular,  the  truncated  infinitesimal  generator  matrix  is 


Q 


T 

n 


-r(0)  1(1) 

r(0)  — r(l)  —  1(1)  1(2) 

r(l) 

l(n  —  1) 

r(n  —  2)  —  r(n  —  1)  —  l(n  —  1) 


and  has  the  following  explicit  LU  factorization 
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'  -1 

'  Ho)  -i(i) 

L  = 

1 

U  = 

Hi) 

-1 

‘  * .  —  l(n  —  1) 

1  _1  . 

r(n  -  1) 

To  find  the  stationary  distribution  we  solve,  in  succession,  the  triangular 
systems 


Lzn  —  en 
U^n  — 

Since  L  is  unit  lower  triangular  we  see  that  zn  =  en  (indeed,  in  general 
one  need  only  know  U)  and  in  is  just  the  solution  of 

Ufn  —  en 

Backward  substitution  yields 


r  f(l)f(2)-/(n)  /(2)f(3).../(n) 

.r(0)r(l)...r(n)’  r(l)r(2)...r(n)’ 

This  implies  that 


l(n)  1 

r(n  —  l)r(n)’  r(n). 


n+1  IN) 

p(k)  =  p(n+  1)  n  rf  •_-!) 

j=*+l  U  7 

Setting  k  =  0  and  rearranging  yields 


(7) 


p(n+l)  =  p(0)n 

which  is  a  well-known  formula  for  the  stationary  distribution. 

Note  that  this  formula  does  not  give  the  values  of  the  stationary  distri¬ 
bution  since  it  involves  the  unknown  scaling  factor  p(0).  However,  it  does 
allow  us  to  determine,  exactly,  the  shape  of  the  stationary  distribution.  We 
can  make  a  common  approximation  to  the  stationary  distribution  in  the 
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following  way.  Assume  that  5  is  the  probability  that  the  population  size  is 
less  than  n.  Then  eTpn  =  S  and  hence 

5  f 

Pn  ~  eTfnn 

In  practice,  one  takes  n  to  be  large  so  that  S  is  close  to  1.  Then 

1  f 

Pn  ~  eTfn  *n 

4  Populations  with  Multiple  Births 

Now  consider  populations  in  which  multiple  births  can  occur.  Although  one 
approach  to  problems  of  this  type  is  to  re-scale  the  birth  rate,  r(n),  in  some 
appropriate  way  and  solve  the  single  step  problem,  it  is  not  difficult  to  derive 
the  solution  using  traditional  methods.  We  show  that,  as  before,  the  matrix 
approach  gives  the  formal  solution  when  solved  analytically. 

Let  r(n,  k)  be  the  rate  at  which  births  of  k  individuals  occur  given  that 
the  population  size  is  n.  Equation  2  becomes 

(oo  \  n 

^2  r(n,  k )  +  l(n)  J  p(n)+J^  r(n-k,  k)p(n-k)+l(n+l)p(n+l)  (8) 
k= l  /  k=\ 

The  infinitesimal  generator  is  an  infinite  lower  Hessenberg  matrix  whose 
elements  are  given  by 

-(££i  r(i,k)  +  l(i))  if  i=j 

r(j,i-j)  if  i>3  /o', 

l(j)  if  i  +  1  =  j  y  ’ 

0  otherwise 

The  truncated  system  has  the  same  form  as  before  except  that  QT  is 
now  lower  Hessenberg  instead  of  tridiagonal.  In  particular,  letting  pn  = 
l(n)p(n)%l  we  have 

Qn^n  =  -e„  (10) 

As  before,  we  only  need  to  know  U  from  the  LU  factorization  of  Qr 
to  find  f„.  In  block  form,  the  factorization  of  the  truncated  infinitesimal 
generator  matrix,  an  n  x  n  lower  Hessenberg  matrix,  looks  like 
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1  0T  1  [  pej  }  _  r  a  fief  ' 

-r  L  0  U  ~  r  QT 

la  JL  JL 

Where  LU  is  the  LU  factorization  of  the  Schur  complement  which  is  an 
n-lxn-1  lower  Hessenberg  matrix.  In  particular 

LU  —  QT  -  —ref 
a 

Note  that  U  is  an  upper  bidiagonal  matrix  whose  first  super- diagonal 
has  elements  Ut<t+i  =  l(i).  Let  <7,-  =  UlA  be  the  diagonal  elements,  these  can 
be  found  sequentially  in  the  following  way. 
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And  setting  k  =  0  and  rearranging  yields 

n+1  <7 

p{n  +  1)  =  p(0)  I] 

5  Processes  with  Multiple  Births  and  Multiple 
Deaths 

In  the  cases  considered  so  far,  the  infinitesimal  generator  matrix  is  lower 
Hessenberg  and  the  solution  algorithms  are  equivalent  with  the  classical 
solution  by  recursion  algorithm  for  finding  truncated  solutions.  The  deriva¬ 
tions  we  have  shown  are  different  than  the  classical  approach  but  are  useful 
because  they  will  allow  us  to  look  at  more  general  birth  and  death  models 
in  a  natural  way.  We  now  consider  a  process  in  which  both  multiple  births 
and  multiple  deaths  are  allowed.  We  will  assume  that  jumps  as  large  as 
±A'  occur  in  both  directions  (it  is  straightforward  to  extend  what  follows  to 
processes  where  the  maximum  possible  number  of  births  is  not  the  same  as 
the  maximum  possible  number  of  deaths).  The  stochastic  balance  equation 
is 


0  =  J2  (r(n»  fc)  +  Kni  k))p{n)  +  r(n  ~  k,  k)p(n  —  k)  +  l(n  +  k,  k)p(n  +  k)} 
k= i  ,  „ 

(11) 

The  infinitesimal  generator  matrix  is  banded,  with  half-bandwidth  R , 
and  of  the  form 


-ZLi  (r(i,k)  +  l(i,k)) 
r(j,  i  ~  j ) 
KhJ-i) 

0 


if  i  =  j 
if  i>  j 
if  i  <  j 

otherwise 


The  truncated  system  can  be  written 


(12) 


QnPn  +  =  0 

where  p(n)  =  [  p(n)  p(n  +1)  ...  p(n  +  K  —  1) 

elements 


T  and  €  $lnxK  has 
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s(n)  _  |  l(j,  j+n-i)  if  i  >  n  -  K  +  j  -  1 
0  otherwise 

We  shall  call  the  homogeneous  complement  of  Q £  in  QT . 

Rearranging  the  truncated  system  yields 

QTnPn  =  -S<">  p<"> 

Now,  let  F  be  the  solution  to 

QTnF  =  -S™ 

Then  pn  =  Fp(n)  and  we  see  that  pn  is  in  the  range  of  F.  If  F  is  rank 
one  then  we  can  get  bpn  up  to  an  unknown  scaling.  This  is  the  essence  of 
solution  by  recursion  since  in  those  cases  F  surely  has  rank  one. 

Based  on  the  preceding  analysis,  we  propose  the  following  algorithm  to 
find  p(0),p(l),...p(iVo)  for  a  birth  and  death  process  that  can  have  from  1 
to  K  deaths  in  each  time  step.  We  assume  that  Nq  >  K . 

1.  Let  n  —  N0. 

2.  Using  the  truncated  infinitesimal  generator  Q ^  and  its  homogeneous 
complement  solve  Q^Fn  —  —S^. 

3.  Compute  the  singular  value  decomposition  of  Fn,  that  is  UT,Vt  =  Fn 

4.  Construct  the  approximation  pn  «  aui  where  a  is  some  unknown 
constant  and  ui  is  the  first  column  of  U, 

5.  If  cr[n\  the  largest  singular  value  of  Fn  is  sufficiently  greater  than 
then  stop  and  accept  the  approximation.  Otherwise  set  n  :=  n  + 1  and 
return  to  step  2. 

It  is  possible  to  update  Fn  and  it’s  singular  value  decomposition  quickly. 

6  A  More  General  Formulation 

The  method  that  has  been  developed  above  can  be  put  into  a  much  more 
general  framework.  We  begin  by  noting  that  any  Markov  process  whose 
infinitesimal  generator  is  banded  is  a  Quasi-Birth-Death  (QBD)  process. 
Provided  we  choose  the  block  size  correctly  (it  must  be  no  less  than  the 
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bandwidth)  the  infinitesimal  generator  can  be  written  as  a  block  tridiagonal 
matrix 

Do  A\ 

B\  D\  A2 

1 T 

If  we  partition  p  =  |  7T0  tti  ...  so  that  it  is  compatible  for  block 
multiplication,  then  the  stochastic  balance  equations  can  be  written  in  the 
following  form 

£>0*0  +  ^4.1 7T 1  =  0 

and 

B{  i — 1  "I”  D  i  7T{  4"  ^ 

for  i  =  1,2, ... 

There  are  two  common  approaches  to  problems  of  this  type.  First,  if  the 
Markov  process  represented  by  this  matrix  is  nearly  completely  decompos¬ 
able  (NCD),  that  is,  if  the  off-diagonal  blocks  are  sufficiently  small  in  some 
sense,  one  can  disregard  them  and  assume  the  Markov  chain  is  completely 
decomposable  (some  justification  for  this  can  be  found  in  [3]).  Then  7r are 
solutions  to  D{-K{  =  0  and  can  be  solved  for  individually.  This  yields  ap¬ 
proximations  to  the  segments  of  p  which  must  then  be  carefully  assembled 
to  get  the  stationary  distribution  (see  [4]  for  a  beautiful  treatment  of  these 
methods).  The  problem  with  this  approach  is  that  it  is  hard  to  tell  how  good 
the  approximations  to  the  t,  are  since  it  is  not  clear  what  effect  disregarding 
the  A{  and  B%  will  have  on  the  solution. 

A  second  approach  is  stochastic  complementation  ([2]).  This  method 
also  computes  the  segments  /T,  but  does  so  exactly  using  block  Gaussian 
elimination.  This  method  does  not  throw  anything  away  so  it  is  exact  (at 
least  on  paper).  Unfortunately,  the  method  is  quite  costly. 

Our  method  can  be  extended  quite  naturally  to  this  more  general  frame¬ 
work.  In  particular,  given  the  block  tridiagonal  QT  shown  above,  we  define 
the  homogeneous  complement  of  £>;  in  QT  to  be  S\  =  j  1  A ;  ,  where 

is  a  matrix  composed  of  only  the  non-zero  columns  of  A,,  B,  is  defined 
similarly,  and  Bo  is  taken  to  be  a  zero  matrix.  We  then  solve 
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DiFz  —  Si 

and  approximate  i r;  with  the  left  singular  vector  associated  with  the  largest 
singular  value  of  F{.  This  is  a  simple  process  and  allows  us  to  estimate  of 
the  quality  of  our  approximations  by  examining  the  ratios  of  the  largest 
and  second  largest  singular  values  of  each  F{.  This  method  is  embarassingly 
parallel. 

We  can  use  this  method  to  generate  a  starting  guess  for  those  algorithms 
known  generally  as  Iterative  Aggregation/ Disaggregation  (IAD).  These  are 
efficient  multi-grid  like  methods  and  include  the  well-known  KMS  [1]  and 
Takahashi  [5]  algorithms.  We  can  also  develop  an  adaptive  variation  by 
using  the  algorithm  of  section  5  to  solve  for  7To  first.  We  then  apply  this 
same  algorithm  to  solve  for  each  succesive  segment  using  the  more  general 
definition  of  the  homogeneous  complement.  This  allows  us  to  vary  the  block 
sizes  adaptively  so  that  each  estimate  of  a  segment  will  be  a  good  one. 
This  approach  lacks  parallelism  but  would  give  both  an  initial  guess  and  a 
partitioning  for  IAD  algorithms. 

7  An  Example 

As  an  example  we  look  at  a  simple  birth- and- death  process  in  which  as  many 
as  two  births  or  deaths  can  occur  in  the  fundamental  time  unit.  The  specific 
model  we  will  look  at  is  characterized  by  the  following  transition  rates 


r(n,  1)  =  0.237(n+  i)e~°*0165n 

r(n,  2)  =  0.105(71 +l)e“’°-0231n 

/(n,  1)  =  0.088n 

/(7i,2)  =  0.018n 

Notice  that  there  is  a  positive  birth  rate  even  when  the  population  size 
is  zero.  This  migration  term  is  necessary  so  that  0  is  not  an  absorbing  state 
which  would  preclude  the  existence  of  a  stationary  probability  vector. 

Analytically,  we  need  to  find  the  null- vector  of  the  full  infinitesimal  gen¬ 
erator  Q,  an  infinite  matrix.  In  practice  it  is  sufficient  to  find  the  eigen¬ 
vector  associated  with  the  smallest  magnitude  eigenvalue  of  Qn.  As  n  gets 
large,  this  eigenvector  should  converge  to  the  stationary  distribution,  after 
appropriate  normalization.  For  this  example,  when  n  =  200  the  smallest 
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eigenvalue  has  magnitude  roughly  10-13  and  we  see  good  convergence  of  the 
eigenvector.  Of  course,  we  had  to  solve  a  large  eigenvalue  problem  to  get  this 
approximation.  Below  is  a  plot  of  the  first  100  elements  of  the  eigenvector 
(past  100  the  distribution  dies  out) 


We  now  explore  the  application  of  the  methods  described  in  this  paper  to 
this  example.  First  of  all,  working  with  the  matrix  truncated  to  the  first  100 
equations  we  take  the  homogeneous  complement,  solve  and  return  the  left 
singular  vector  associated  with  the  largest  magnitude  singular  value.  The 
ratio  of  the  largest  to  the  second  largest  singular  value  is  roughly  3  X 105.  The 
approximation  is  almost  identical  to  that  found  by  solving  a  full  200  x  200 
eigenproblem.  Below  is  a  plot  of  the  difference  between  the  approximations 
after  appropriate  normalization. 
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Notice  that  the  difference  never  exceeds  2.5  X  10-7  in  magnitude.  More¬ 
over,  the  relative  error  in  the  approximation  for  all  n  such  that  p(n)  >  .005 
never  exceeds  3  X  10-12.  Note  also  that  we  did  not  have  to  solve  a  200  X  200 
eigenproblem  to  get  this  solution.  We  only  had  to  solve  a  single  linear  sys¬ 
tem,  with  a  100  x  100  matrix  and  a  100  x  2  right  hand  side,  and  compute 
the  singular  value  decomposition  of  a  100  X  2  matrix. 

If  we  truncate  the  infinitesimal  generator  so  that  it  covers  population 
sizes  ranging  from  40  to  90  only  we  find  that  the  ratio  of  the  two  largest 
singular  values  of  F  is  roughly  36.4.  After  appropriate  normalization,  the 
difference  between  this  segment  and  the  full  solution  is  always  less  than 
6  X  10-4,  we  see  this  plotted  below. 
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And  finally  we  show  a  plot  of  the  computed  segment  (dotted  line)  with 
the  segment  from  the  full  solution. 


We  see  that  the  method  works  quite  well  in  this  case  and  is  much  less 
costly  since  the  SYD  computation  involves  a  very  small  matrix. 
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